import numpy as np


def steepest_descent(A, b, x0, tol=1e-4, max_iter=1000):
    n = len(b)
    x = np.zeros(n) if x0 is None else x0.copy()
    history = []
    for k in range(max_iter):
        g = A @ x - b  # 梯度 ∇f(x) = Ax - b
        d = -g
        alpha = (d @ d) / (d @ A @ d)  # 最优步长
        history.append((x.copy(), d, g, alpha))
        x = x + alpha * d

        if np.linalg.norm(g) < tol:
            break
    return x, history


if __name__ == "__main__":
    A = np.array([[2.0, 0.0, 0.0], [0.0, 1, 0.0], [0.0, 0.0, 1]])
    b = np.array([0.0, 0.0, 0.0])
    x0 = np.array([1.0, 1.0, 1.0])
    x_opt, history = steepest_descent(A, b, x0)
    for i in range(len(history)):
        print(
            f"Iteration {i}: x = {history[i][0]}, d ={history[i][1]}, g = {history[i][2]}, alpha = {history[i][3]}"
        )
    print(x_opt)
